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FOREWORD 


This is a progress report on the research project "Numerical Solutions 
of Three-dimensional Navier-Stokes Equations for Closed-Bluff Bodies." The 
period of performance on this research was from January 1 through August 
31, 1986. The project is supported by the NASA/Langley Research Center, 
Hampton, Virginia, and monitored by Dr. Robert E. Smith, Jr., of the 
Analysis and Computations Division (Computer Applications Branch), 
NASA/Langley Research Center. 
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ABSTRACT 

Methods of grid adaption are reviewed and a method is developed with the 
capability of adaption to several flow variables. This method is based on a 
variational approach and is an algebraic method which does not require the 
solution of partial differential equations. Also the method has been 
formulated in a such way that there is no need for any matrix inversion. The 
method is used in conjunction with the calculation of hypersonic flow over a 
blunt nose. The equations of motion are the compressible Navier-Stokes 
equations where all viscous terms are retained. They are solved by the 
MacCormack time-splitting method and a movie has been produced which shows 
simultaneously the transient behavior of the solution and the grid adaption. 
The results are compared with the experimental and other numerical results. 


Graduate Research Assistant, Department of Mechanical Engineering and 
Mechanics, Old Dominion University, Norfolk, VA 23508, AIAA Student 
Member. 

2 Eminent Professor, Department of Mechanical Engineering and Mechanics, Old 
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INTRODUCTION 


In hypersonic flows about bluff bodies, the temperature, pressure and 

density increase almost explosively across the shock wave and the curved shock 
wave is closer to the body. Numerical simulations of this phenomena is a 
great challenge to the computational fluid dynamics researchers. Presently, 
there is a great deal of interest in improving the quality of numerical 
simulation techniques and adaptive gridding is one way to achieve this goal. 

Grid generation is the very first step in the numerical solutions of 
partial differential equations for complex geometries. Basically, grid 
generation is the numerical generation of boundary-fitted curvilinear 
coordinates. The second step is the construction of difference equations for 
partial differential equations. It is apparent that the accuracy of the 

finite difference solution depends on the fineness of the mesh. Therefore, 
the finer the grid, the more accurate the numerical solution will be. Also, 
the accuracy of the solution depends on the resolution of the solution 
gradient. Presence of large gradients causes the error to be large in the 

difference approximation of the derivatives. In the presence of shockwave, 
more artificial diffusion must be added to retain adequate smoothness of the 
solution. Therefore, there is a need for schemes that can resolve these large 
gradients without adding any additional grid points. An adaptive scheme moves 
the grid points to regions of high gradients when locations of these gradients 
are not known as a priori. Also, an adaptive method reduces the total number 
of grid points required to achieve a given accuracy, but it takes more 
computer time. In some instances, the computer time requirement makes this 
method undesirable. The ideas used in the construction of adaptive grid 
techniques are limited by one's imaginations; any scheme that works in the 

sense of providing a better solution is a good one. The ultimate answer to 
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numerical solutions of partial differential equations may well be dynamically 
adaptive grids rather than more elaborate difference representations and 
solution methods [1]. 

Adaptive methods have been used in the solutions of ordinary differential 
equations. Variable-step initial-value problems are solved by adjusting the 
step size as the integration advances in order to control the local truncation 
error [2]. Adaptive methods have been implemented also for solving the 
equations of motion in conjunction with the method of lines [3]. In this 
case, the time step was adjusted automatically to control the local error. 
Similarly, adaptive methods have been to solve boundary value problems [4- 
9], An optimal grid for a two-point boundary-value problem can be determined 
either implicitly or explicitly. In the implicit approach, the weight 

function depends upon the solution. As a result, the original boundary-value 
problem is converted into an augmented system in which the dependent variables 
and the mesh are computed simultaneously. In the explicit approach, the 
weight function does not depend on the solution. Instead, it depends upon a 
previously calculated solution. Even for a linear boundary-value problems, 
the implicit approach requires that one solves a nonlinear two-point boundary- 
value problem. The implicit techniques do not preserve the linear-nonlinear 
character of the original problem. Moreover, even for the nonlinear problem, 
the augmented system is usually more difficult to solve than the original 
problem. On the other hand, the explicit technique preserves the linear- 
nonlinear character of the original two-point boundary-value problem. 

Adaptive schemes are divided into two basic categories: differential and 

algebraic. The differential method is based on the variational approach. 
Brackbill and Saltzman [10-13] have developed a technique for constructing 
adaptive grids using a variational approach. In their scheme, a function 
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which contains a measure of grid smoothness, orthogonality and volume 
variation is minimized by using the variational principal. The smoothest grid 
can be generated by solutions of Lapl.ace equations which are better known as 
elliptic systems [14, 15]. This approach ignores the effects of orthogonality 
and requires too much CPU time. The method has been modified for better 
efficiency by dropping the second derivative terms in one coordinate direction 
[16]. This makes the equations parabolic and, therefore, they can be solved 
by marching techniques. A method which considers the orthogonality and volume 
variation is developed by Steger and Sorenson [17], This method is widely 
known as the hyperbolic method; it can be solved by any non-iterative marching 
techniques. The variational approach provides a solid mathematical basis for 
the adaptive method, but the Euler-Lagrange equations must be solved in 
addition to the original governing differential equations. On the other hand, 
the algebraic method requires much less computational effort, but the grid may 
not be smooth. 

Rai and Anderson [18-23] have developed an algebraic technique where the 
grid movement is governed by estimates of the local error in the numerical 
solution. This is achieved by requiring the points in the large error regions 
to attract other points and points in the low error region to repel other 
points. Nakahashi and Deiwert [24] have formulated an algebraic method which 
is based on the variational principal. A spring analogy is used to redistri- 
bute the grid points in an optimal sense to reduce the over all solution 
error. In this case, operator splitting and one sided controls for the 
orthogonality and smoothness are used to make the method practical, robust and 
efficient. Dwyer [25-28] has used an adaptive method in which the points are 
moved along one set of the original coordinate lines in response to the 
evolving gradients in the physical solution. The analysis showed that the 
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percentage change in a dependent variable can be determined as a priori. An 
order of magnitude improvement in speed was obtained, but some problems with 
excessive skewness were encountered. 

Generally, dynamic adaption can be performed in two ways, one is to keep 
the computational space fixed and include the grid speed in the flow field 
equations. This is an ideal method to use for unsteady flows. The second way 
is to set the grid speed equal to zero and interpolate the solution onto the 
new grid after each adaption. This is equivalent to solving a sequence of 
boundary-value problems, which is an economical way to treat steady flows 
whose solutions are approached asymptotically. It is generally sufficient to 
adapt just a few times during the course of the computation. In this 
approach, the grid distribution at time N+l is determined from time N. Dwyer 
adapted the grids after each integration step or after a selected number of 
steps [25]. However, the grid speed can be obtained by postulating a law 

which is based on some solution properties. Then, these equations can be 
integrated with the governing partial differential equations to yield a new 
grid distribution [29]. The advantage of this technique is that the grid 
location and grid speed are time accurate. 

The literature survey indicates most techniques adapt to just one 
variable. This means that the weight function is based on the solution of one 
variable only. However, the solution of equations of motion produce several 
dependent variables. Viscous-hypersonic flow over a blunt-body have large 
gradient in pressure, velocity, etc in different part of the flow field. 
There is a large gradient in pressure near the shock region; at the same time, 
there is a large gradient in velocity near the solid body. Therefore, there 
is a need for the development of efficient grid adaption method which utilizes 
several variables simuntanously. 
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METHODS OF GRID ADAPTION 

One reason to use grid adaption is to minimize the error over some domain 
by rearranging the grid. Calculus of variation can be used to perform this 
minimization. In general, a weighted integral which is a measure of some grid 
or solution property over some domain can be expressed as 

I = / W d V, (I) 

V 

where W is the weight function and it is minimized. The selection of W may 
vary from problem to problem. There is a collection of definitions for W in 
Ref. 30. The weight function can be based on grid properties such as cell 
volume, the average of the square of diagonal lengths, the cell area/volume 
ratio and/or cell skewness. There exists a differential equation which 

minimizes the integral I in Eq. (1). This differential equation is called 

Euler-Lagrange differential equation [31] and it constitutes the grid 

generation system. The Euler-Lagrange equations can be found in Refs. 10-13. 

Multidimensional Adaption 

Brackbill and Saltzman [10-13] have developed a technique based on a 
variational approach. In their scheme, a function which contains measures of 
grid smoothness, orthogonality, and volume variation is minimized. To 

maximize the smoothness of the grid, the following integral must be minimized 

3 

I = / 2 V S 1 • V? 1 dV. (2) 

s v i=i 

This is simply the sum of the squares of cell-edge lengths. Similarly, the 
orthogonality can be acquired by minimizing the integral I Q given by 

I « J g 3/2 (VC 1 ‘VC 1 ) dV. 

0 V 


(i,j,k) cyclic 


(3) 
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This integral vanishes for an orthogonal grid. The inclusion of the Jacobin 
of the transformation in the weight function is somewhat arbitrary and causes 
orthogonality to be better in the regions with the larger cell. The concen- 
tration or cell volume variation can be obtained by minimizing the integral 

I = / W J dV, (4) 

w v 

where W is a specified weight function. This causes the cells to be small 
where the weight function is large. The grid generation system which provides 
smoothness, orthogonality and concentration is obtained by minimizing the 
total integral I which is a linear combinations of I s ,I 0 and I w 


I 


I + X 
s o 


I + X I 


w w 


(5) 


The competing features such as smoothness, orthogonality, and cell volume 
variation can be stressed by proper choice of the coefficients X Q and X w> 
For example, a large X q will result in a nearly orthogonal grid at the cost 
of the smoothness and concentration. The Euler-Lagrange equations for the 
sums of those individual integral form the system of partial differential 
equations from which the coordinate system is generated. These equations are 
quasi-linear, second-order partial differential equations with the coeffi- 
cients which are quadratic functions of the first derivatives [12]. This 
variational formulation is equivalent to Winslow's method [14] where \ Q and 
X are set equal to zeros. The Euler equations are those given by Winslow, 
and their solution maximizes the smoothness. This is also used by Thompson et 
al. [15], The additional terms alter other characteristics of the mapping in 
a similar way. The cell-size variation and skewness can be controlled by 
proper selection of I $ , I Q and I w . The use of a variational approach provides 
a solid mathematical basis for the adaptive gridding. But, the Euler-Lagrange 
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equations must be solved in addition to the governing equations of fluid 
motion. 

One Dimensional Adaption 

The Euler-Lagrange equation for Eq. (5) is a general and is capable of 
adapting grids in multidimensions simultaneously. When the solution varies 
predominately in a single direction, one-dimensional adaption can be applied 
with the grid points constrained to move along one family of fixed curvilinear 
coordinate lines. The fixed family of lines is established by generating a 
full multidimensional grid using the grid generating techniques. The points 
generated for these initial grids, together with some interpolation procedure, 
e.g., cubic or linear interpolation, serve to define the fixed lines along 
which the points will move during the adaption. This is done explicitly; 
therefore there is no need to solve any differential equations. 

A technique called equidistribution is developed to improve the solutions 
of boundary value problems [4-9], This has proven to be effective and 

efficient. This technique is used to minimize the error by redistributing the 
grid such that a weight function is constant over each interval. The Euler- 
Lagrange equation is 


Xjj W = constant. (6) 

This minimizes the following integral 

1 ? 

I, = / WU) K d? . (7) 

1 0 * 

This equation represents the energy of a system of springs with the spring 
constant W(£), spanning each grid interval. The weight function is 
associated with the grid points themselves, not with their locations. An 


8 


alternative viewpoint results by integrating over x instead of £, i.e. 
summing over the grid intervals rather than over the grid points. This can be 
expressed as 

1 4 

l 2 = f Q TTTxT dx ' (8) 

The Euler-Lagrange equation for this is given by Eq. (6). The quantity C x is 
considered to represent the point density. This variational problem 
represents a minimization over the density of the grid points subjected to a 
weight function. This can produce smooth point distribution. Here the weight 
function W(x) is associated with the location of points. If the weight 
function is associated with the points themselves rather than their locations, 
then W = W(5). The integral for which Eq. (6) is the Euler-Lagrange 
equation for the following integral 

*3 = ^ dx> (9) 

£ x is a measure of the smoothness of the point distribution, with the 
emphasis placed on the smoothness in the certain regions. This is inversely 
proportional to the weight function W(£). Equation (6) is the Euler- 

Lagrange equation for the integrals in Eqs. (7-9) which can be written as 

X 

5(x) = / W(5) dx. (10) 

0 

Equation (10) can be written in terms of the arc length as 

s 

S(s) = / WU) dx. (11) 

0 

The weight function is used to reduce the point spacing where W is large. The 
weight function should be set to some measure of the error. White [5] has 
suggested the following form of the weight function 
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r 

W = [a + |u {n) | 2 ] 1/2n , 

where a is constant. With n=l and a=0, this becomes 


l 

A combination of Eqs. (6) and (13) yields 
| = Constant. 


( 12 ) 


(13) 


(14) 


This choice replaces the points so that the same change in the solution occurs 
at each grid interval. This is simply the solution gradient. Taking n=l and 
a=l yields 



(15a) 


A combination of Eqs. (6) (15a) results in 



Constant. 


(15b) 


This produces a uniform distribution of arc lengths on the solution curve. 
White's results [5] indicated that the arc length form was favored. The 
disadvantages of this method is that the weight function near the solution 
extreme, i.e. U x =0 locally, are treated as a flat region. Concentration near 
the solution extreme can be achieved by incorporating some effect of the 

second derivative (U xx ) into the weight function as [25] 

W = 1 + a f(U x ) p g(U xx ), (16a) 


where a and {3 are positive parameters. Therefore, Eq. (10) can be written as; 
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x 

/ W(x) dx 

5 (x) = -j . (16b) 

/ W(x) dx 
0 

With the second derivative terms included, the value of a, p must be 
continually updated to keep the same relational emphasis or concentration 
according to this form. Therefore, a system of 2x2 should be solved for each 
fixed grid line [27-28]. It will be shown later that through the reformula- 
tion of Eq. (16b), the matrix inversion can be entirely avoided. 

One-Dimensional Adaption With Several Variables 

The solution of the flow equations consists of several variables. 
Therefore, the weight function should be a function of more than one 
variable. It is desirable to devise a scheme in which grids can adapt to 

several variables with control on the magnitude of adaption for each 
variable. In the case of high speed flow, velocity has large gradients in 

some regions whereas pressure is constant or vice versa. In general the 
weight function can be expressed as 


N 

W = 1 + Z b. f. (17) 

i=l 1 1 


where N is the number of variables, b^ is constant, fj is some variable or its 
derivative, and 1 is for uniformity. A substitution of Eq. (17) into Eq. (11) 
results in 


5(S) 


N 

S + Z b. F. (S) 
i=l 1 1 

R 


s + z b. f. (s ; 

max l l max 


(18) 
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where 


F i 


S 


(S) = / 

0 


fj (t) dt . 


It should be noted that Eq. (18) is for adapting along a fixed grid line. The 
quantities and fj should be positive to ensure the monotoncity of £(s). 
In order to keep the same relative emphasis on the concentration along each 
grid line, b.,- should be computed based on some percentage of the grid being 
allocated to each variable. The percentage of grid points assigned to a 
particular function fj can be expressed as 


R. = 
J 


“j VW 


N 

S + E b. F. (S ) 
max i i max 


(19) 


j=l,2. 


N. 


A rearrangement of Eq. (19) results in 


[ V { V ■ V • 


where 


A.. - R . Ilife*! 

J1 J TO 


1 * J 


( 20 ) 



i = J, 


C j = " R J S max /F j (S max } * 


Therefore, along each fixed grid line a matrix of NxN (N is number of 
variables) should be inverted. This can be avoided by the reformulation of 
Eq. (18). The crucial steps are outlined here. Equation (18) can be 


rewritten as 
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5 (S) = 


R 

S + Z b. F. (S ) 

max j=i 1 1 max 


N 

+ Z 
i=l 


b. F. (S) 


N 

S + Z b . F . (S ) 
max . , j j max' 


( 21 ) 


Similarly, Eq. (19) can be rewritten as 


Vnrj^V^tVWl' (22 > 

J FTia a 1 — 1 


A substitution of Eq. (22) into (21) results in 


N R_. 


5(S) = 


. 5 * [s -« + ji. 


(S) 


+ [$ mav + 2 b. F. (S )] Z R. 
L max 1 l max J j 


max 


(23) 


A summation of Eq. (19) over all j yields 


N 

Z R. = 

j-1 J 


Z b. F. (S ) 
j j max' 


j-1 


N 

S + Z b. F. (S ) 
max -j=i 1 1 max 


Rearrangement of Eq. (24) results in 


N 

Z b. F. (S ) = 
l i max 


N 

S max 1 f 1 R i 

R 

1 - Z R. 
1=1 1 


A substitution of Eq. (25) into Eq. (23) yields 


(24) 


(25) 


F, (S) 


£(S) = [l 2 R-] + Z R. j — j-* 

max i=l i=l h i max' 


( 26 ) 
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This reformulation avoids the need for continuous updating the b^'s to keep 
the same relative emphasis on concentration. Applying this equation, grids 
can be adapted with more than one variable without any need for matrix 
inversions. 

GOVERNING EQUATIONS OF MOTION 

The governing equations for a thermal fluid system are the conservation 
of mass, momentum, and energy. These equations are developed for an arbitrary 
region assuming the system is in continuum. Equations of motion for a viscous, 
compressible, unsteady, heat conducting fluid can be written as [32] 

Continuity: + V • (p u) =0 , 

Momentum: H • (p y u ■ t) = 0 , (27) 

— — — = 

Energy: ■^+V»(Eu + q- u : c)=0. 

where E is the total energy per unit volume given by (E=e+u»u/2) and e is 
the internal energy per unit volume. The fluid is assumed to be Newtonian and 
the bulk viscosity is neglected. The viscosity is computed from the 
Sutherland law of viscosity. For simplicity, equations of motion can be 
written into a compact vector form as 
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For the sake of generality, governing equations are transformed from a 
physical domain into a computational domain. Resulting equations are: 


dll 

■5t 


+ 



r dG SH 
* -Sri 


) + 



r 5G 5H 


) - o. 


(29) 


The transformation coefficient can be computed from a functional relation 
between the computational coordinates and the physical coordinates. 

METHOD OF SOLUTION 

A time marching method is used to compute the solution. This allows us 
to capture the possible transient features. The method employs an explicit 
second-order accurate time-split predictor-corrector algorithm [33]. In a 
compact form, it can be expressed as 


n+1 

U j>k = [L^ (Ay] [L c (At c ) ] ly (Ay] , (30) 

where 

The method' has a time-step stability limit, but there is no rigorous nonlinear 
stability analysis available for it. However, there is a conservative time- 
step which is based on a linearized form of the equations. It can be 
expressed as 

At<Mn [M.fl.c [1 ♦ * ]] 1/2 (31) 

L Ail AC L Aq? 
where c is the local speed of sound. 

In the hypersonic region, there exits a large gradient which requires a 
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very fine mesh to resolve it. Most center difference methods admit a solution 
which has sawtooth or plus-minus waves with the shortest wave length that the 
mesh can support. In the of case nonlinear problems, these short waves 
interact, vanish, and reappear again as distorted long wave or oscillation. 
These oscillation eventually blow up the solution, if they are not resolved. 
The oscillations of "low frequency" can be suppressed by adding a fourth order 
damping. A common damping used is the pressure dampening. This can be 
expressed in computational coordinates as 


a+ c3 d pi 
a JL At JL 6 Jl "3S 


+ c a p 


au 


~W 




where is contra variant velocity. 


(32) 


RESULTS AND DISCUSSIONS 

The main objective of this study is to investigate finite difference 
methods in which the mesh network adapts to the solution dynamically to obtain 
an accurate solution for hypersonic flow. A computer program has been written 
to utilize Eq. (26) for grid adaption. Presently, this code is being run on 

the Network Operating System (NOS) and the Vector Processing System (VPS) at 

NASA Langley Research Center. Hypersonic flow over a blunt nose is a typical 
test case in computational fluid dynamics. This problem inherits a detached 
shock which should be accurately resolved. The location and the magnitude of 
the shock are not known as a priori; therefore the grid should be adapted as 
the solution progresses. This problem is used to analyze and verify the 
adaptive method. Equations of motion are solved by the MacCormack method [33] 
for a hypersonic flow over a small-radius blunt-body with an inclined-plate 
after the body (Fig. 1). The blunt leading edge is a part of the panel holder 

which was tested at Langley Research Center [34]. Results are obtained at 
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following conditions: M m =6.8, Re oo =220,000, P^-9.26 lb/ft 2 , 11^=6510 

ft/ sec, y=1.38, R=1771 ft 2 /sec 2 /°R and wall temperature of 540°. 

Two tests have been performed: static adaption and dynamic adaption. For 
static adaption, the solution has been obtained by fixed grids shown in Fig. 
2. The grids are adapted to two variables, the first and the second 
derivatives of pressure, and the results are shown in Fig. 3. In this case, 
twenty percent of the grids are allocated to first and second derivatives of 
pressure (R1=R2=20%). For the same case. Fig. 4 shows adaption with 
R1=R2=502. It is noted that in Fig. 4 entire grids are allocated to the first 
and the second derivatives of the pressure. This is explains the large voids 
in the constant pressure regions. Figures 3 and 4 lack grid resolution in the 
vicinity of the solid boundaries. This is due to the constant pressure near 
the solid boundaries. But Eq. (26) can adapt to several variables. Figure 5 
shows the grids which are adapted to two variables, pressure and velocity. 
The weight function consists of the first derivatives of pressure and 
velocity and the second derivative of pressure. Twenty percent of grids are 
allocated to each function and forty percent of grids are for the uniformity. 
This avoids the creation of any void. It is noted that grids are clustered 
near the shock and the solid body. Therefore, it is possible to resolve the 
pressure as well as the velocity gradients in the boundary layer region. In 
chemically reacting flows, some of the grids can be allocated for resolving 
the gradients of the chemical species. 

The procedure has been applied dynamically to the same problem. Figure 
6a shows the initial grid distribution for this problem. Figures 6 and 7 show 
sequences of grid distributions at different times. In this case, grids are 
adapted to six variables; these are the first and second derivatives of 
pressure, temperature and velocity. Ten percent of grids are allocated to 
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first derivatives equally and five percent to their second derivatives. 
Fifty-five percent of grids are also allocated to the uniformity. A movie has 
been produced of this work which shows the dynamic adaption; a few frames are 
shown in Figs. 7a-7b. It demonstrates how grids are attracted toward high 
gradient region and repel from low gradient regions. The full paper will 
cover formulations, boundary conditions, grid generation and the vectorizing 
of the proposed scheme in more details. Also, the timing of the scheme will 
be reported. Currently, this work is being extended to the blunt nose with the 
after body. 
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Fig. 6 Adapted Grid (Dynamic) 
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Fig. 7 Adapted Grid (dynamic) 
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